function est_m = estimator_s_A_p4s3(est, estinfo)
% fixed effects
% cubic in size

if ~isfield(estinfo, 'FitMethod');
    estinfo.FitMethod = 'ML';
end

nobs = size(est.YY,1);
ncity = length(est.xtval);
nyear = 6;

temp_y = est.YY;
% temp_x = [est.XX, est.logDC];
% temp_yd = [year_dummy2, year_dummy3, year_dummy4, year_dummy5, year_dummy6];
temp_yd = est.YD;

temp_XX = est.XX;
temp_ZZ = est.ZZ;
temp_ZZ2 = est.ZZ2;
temp_logDC = est.logDC;


islog = estinfo.islog;
if islog == 0
    temp_d = [est.DD];
elseif islog == 1
    temp_d = [log(est.DD)];
elseif islog == 2
    temp_d = [log(1+est.DD)];
end

xtflag = est.xtflag;
xtval = est.xtval;

xtflag2 = est.xtflag2;
xtval2 = est.xtval2;

xtflag3 = est.xtflag3;
xtval3 = est.xtval3;

ZZ = est.ZZ;
ZZ2 = est.ZZ2;


%% Model 
Options = [];
Options.Display = 'iter';
Options.TolFun = 1e-12;
% Options = optimoptions('fminunc', 'TolFun', 1e-15, 'TolX', 1e-15, 'MaxFunEvals', 2000);

temp_x = [temp_XX,temp_yd, temp_ZZ(:,2), ones(nobs,1), temp_d, temp_logDC, temp_d.^2, temp_d.^3, temp_d.^4, ...
    (temp_d).*ZZ(:,2), (temp_d.^2).*ZZ(:,2), (temp_d.^3).*ZZ(:,2), (temp_d.^4).*ZZ(:,2), ...
    temp_XX(:,14).^2, temp_XX(:,14).^3 ...
    ];
temp_m = fitlmematrix(temp_x,temp_y,{[ones(nobs,1), temp_d, temp_d.^2, temp_d.^3, temp_d.^4], temp_yd},{xtflag,xtflag},...
    'CovariancePattern',{'FullCholesky','Diagonal'},...
    'RandomEffectPredictors',{{'ej1','ej2','ej3','ej4','ej5'},{'t2','t3','t4','t5','t6'}},'RandomEffectGroups',{'city','time'},...
    'OptimizerOptions',Options, ...
    'FitMethod', estinfo.FitMethod, ...
    'Verbose', 0);
%     'Optimizer', 'fminunc', 'Verbose', 1);
clear est_m
est_m.model = temp_m;

[bet_m] = fixedEffects(temp_m);
est_m.bet = bet_m(1:14);  %coef
est_m.td  = bet_m(15:19); %time-dummy
est_m.a0  = bet_m(21);    % intercept
est_m.a1  = bet_m(20);    % area on intercept
est_m.d10 = bet_m(22);    % distance from the center, linear term
est_m.d20 = bet_m(24);    % distance from the center, quad term
est_m.d30 = bet_m(25);    % distance from the center, cubic term
est_m.d40 = bet_m(26);    % distance from the center, fourth term
est_m.c10 = bet_m(23);    % coastal line, linear term

est_m.d11 = bet_m(27);    % distance from the center, linear term
est_m.d21 = bet_m(28);    % distance from the center, quad term
est_m.d31 = bet_m(29);    % distance from the center, cubic term
est_m.d41 = bet_m(30);    % distance from the center, fourth term

est_m.bet2 = bet_m(31:32); %additional controls for the size effect

[ran_m,ran_m_name] = randomEffects(temp_m);
est_m.ej = ran_m(strcmp(ran_m_name.Group,'city')&strcmp(ran_m_name.Name,'ej1'),:);
est_m.ej2 = ran_m(strcmp(ran_m_name.Group,'city')&strcmp(ran_m_name.Name,'ej2'),:);
est_m.ej3 = ran_m(strcmp(ran_m_name.Group,'city')&strcmp(ran_m_name.Name,'ej3'),:);
est_m.ej4 = ran_m(strcmp(ran_m_name.Group,'city')&strcmp(ran_m_name.Name,'ej4'),:);
est_m.ej5 = ran_m(strcmp(ran_m_name.Group,'city')&strcmp(ran_m_name.Name,'ej5'),:);

est_m.et = reshape(ran_m(strcmp(ran_m_name.Group,'time'),:), nyear-1,[])'; %city by nyear-1
est_m.et = [zeros(ncity,1), est_m.et];

alp_j = est_m.a0 + est_m.a1*ZZ2(:,2) + est_m.ej;
del_j = est_m.d10 + est_m.d11*ZZ2(:,2) + est_m.ej2;
del2_j = est_m.d20 + est_m.d21*ZZ2(:,2) + est_m.ej3;
del3_j = est_m.d30 + est_m.d31*ZZ2(:,2) + est_m.ej4;
del4_j = est_m.d40 + est_m.d41*ZZ2(:,2) + est_m.ej5;

alp_t = repmat([0; est_m.td]',ncity,1) + est_m.et;

alp_jt = [];
alp_jt_mat  = nan*ones(length(xtval), length(xtval2));
alp_jt_mat2 = alp_jt_mat; %smooth version (extrapolate)
for j = 1:1:length(xtval)
    for t=1:1:length(xtval2)
        temp_xtval = xtval(j)+xtval2(t)/10000;
        temp_sel = xtval3 == temp_xtval;
        if sum(temp_sel) ~=0
            
            %temp_alp = alp_j(j) + alp_t(j,t) + est_m.ejt(temp_sel);
            temp_alp = alp_j(j) + alp_t(j,t);
            alp_jt = [alp_jt; temp_alp];
            alp_jt_mat(j,t) = temp_alp;
            alp_jt_mat2(j,t) = temp_alp;
        else
            alp_jt_mat2(j,t) = alp_j(j) + alp_t(j,t);
        end
    end
end

est_m.del_j = del_j;
est_m.del2_j = del2_j;
est_m.del3_j = del3_j;
est_m.del4_j = del4_j;

est_m.alp_j = alp_j;
est_m.alp_t = alp_t;
est_m.alp_jt = alp_jt;
est_m.alp_jt_mat = alp_jt_mat;
est_m.alp_jt_mat2 = alp_jt_mat2;

%% spatial function
est_m.sp = get_sp2(est_m, islog);

%% posterior moments of random effects
% prior calibration
[psi,mse,stats] = covarianceParameters(est_m.model);

big_cov_ad = psi{1}; %cov(alp,del)
big_cov_at = psi{2}; %cov(alp_jt*)
% big_cov_at = blkdiag(0, big_cov_at);

ndel = size(big_cov_ad,1)-1; %number of deltas
nT = 6; % number of years
var_alp = big_cov_ad(1,1); %var(alp)

% cov_aa = var_alp*ones(size(big_cov_at)) + big_cov_at;
cov_aa = blkdiag(var_alp, big_cov_at);
% cov_ad = repmat(big_cov_ad(2:end,1),1,nT);
cov_ad = [big_cov_ad(2:end,1), zeros(ndel,nT-1)];
cov_dd = big_cov_ad(2:end,2:end);

Sig2 = [cov_aa, cov_ad'; cov_ad, cov_dd];
sig2 = mse;

% "Data" for EB procedure: residual form FE part
if ~isfield(est_m, 'a0')
    est_m.a0 = 0;
end
if ~isfield(est_m, 'a1')
    est_m.a1 = 0;
end
if ~isfield(est_m, 'd10')
    est_m.d10 = 0;
end
if ~isfield(est_m, 'd11')
    est_m.d11 = 0;
end
if ~isfield(est_m, 'd20')
    est_m.d20 = 0;
end
if ~isfield(est_m, 'd21')
    est_m.d21 = 0;
end
if ~isfield(est_m, 'd30')
    est_m.d30 = 0;
end
if ~isfield(est_m, 'd31')
    est_m.d31 = 0;
end
if ~isfield(est_m, 'd40')
    est_m.d40 = 0;
end
if ~isfield(est_m, 'd41')
    est_m.d41 = 0;
end
temp_resid = temp_y -temp_XX*est_m.bet -[temp_XX(:,14).^2, temp_XX(:,14).^3]*est_m.bet2...
    - est_m.c10*temp_logDC ...
    ;
%     - temp_yd*est_m.td ... 
%     - (est_m.a0 + est_m.a1*temp_ZZ(:,2)) ...
%     - (est_m.d10 + est_m.d11*temp_ZZ(:,2)) ...
%     - (est_m.d20 + est_m.d21*temp_ZZ(:,2)) ...
%     - (est_m.d30 + est_m.d31*temp_ZZ(:,2)) ...
%     - (est_m.d40 + est_m.d41*temp_ZZ(:,2));


post_y = temp_resid;

post_x = [];
for i = 1:1:ndel
    post_x = [post_x, temp_d.^i];
end
% post_x = [sum(est.YD,2)==0, est.YD, post_x];
post_x = [ones(size(est.YD,1),1), est.YD, post_x];

% Define posterior distribution
V0inv = inv(Sig2);
% posterior moments (bet|sig^2)
postV = @(x,y) (V0inv + sig2^(-1)*(x'*x))^(-1);
% postm = @(x,y)V1*(V0inv*bet0 + sig2^(-1)*(x'*y));
postm = @(x,y,bet0) postV(x,y)*(V0inv*bet0 + sig2^(-1)*(x'*y));


% prior
prior_bet0 = [...
    est_m.a0, est_m.a1; ...
    est_m.td(1), 0;
    est_m.td(2), 0;
    est_m.td(3), 0;
    est_m.td(4), 0;
    est_m.td(5), 0;
    est_m.d10, est_m.d11; ...
    est_m.d20, est_m.d21; ...
    est_m.d30, est_m.d31; ...
    est_m.d40, est_m.d41; ...
];

% RE estimation (BLUP)
nr = size(Sig2,1); %number of random effects
est_m.post.m = zeros(nr,ncity); %posterior mean
est_m.post.V = zeros(nr,nr,ncity); %posterior variance
est_m.post.sig2 = sig2; %variance of measurement error
est_m.post.Sig2 = Sig2;

for j=1:1:ncity
    temp_post_x = post_x(xtflag==xtval(j),:);
    temp_post_y = post_y(xtflag==xtval(j),:);
    %temp_bet0  = zeros(size(Sig2,1),1);
    
	temp_post_z = temp_ZZ(xtflag==xtval(j),:);
    temp_bet0 = prior_bet0*temp_post_z(1,:)';
    
    est_m.post.V(:,:,j) = postV(temp_post_x, temp_post_y);
    est_m.post.m(:,j)   = postm(temp_post_x, temp_post_y, temp_bet0);
    
end


